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Abstract 

In  this  paper,  we  describe  a  method  for  segmenting  fiber 
bundles  from  diffusion-weighted  magnetic  resonance  im¬ 
ages  using  a  locally-constrained  region  based  approach. 
From  a  pre-computed  optimal  path,  the  algorithm  propa¬ 
gates  outward  capturing  only  those  voxels  which  are  locally 
connected  to  the  fiber  bundle.  Rather  than  attempting  to  find 
large  numbers  of  open  curves  or  single  fibers,  which  indi¬ 
vidually  have  questionable  meaning,  this  method  segments 
the  full  fiber  bundle  region.  The  strengths  of  this  approach 
include  its  ease-of-use,  computational  speed,  and  applica¬ 
bility  to  a  wide  range  of  fiber  bundles.  In  this  work,  we  show 
results  for  segmenting  the  cingulum  bundle.  Finally,  we  ex¬ 
plain  how  this  approach  and  extensions  thereto  overcome 
a  major  problem  that  typical  region-based  flows  experience 
when  attempting  to  segment  neural  fiber  bundles. 


1.  Introduction 

Region-based  approaches  to  image  segmentation  consti¬ 
tute  a  key  methodology  for  numerous  applications.  In  these 
approaches,  the  objective  is  to  find  the  segmentation  which 
optimally  separates  features  exterior  to  a  closed  curve  or 
surface  from  features  contained  in  the  interior.  These  ap¬ 
proaches  have  been  shown  to  accurately  segment  datasets 
with  low  signal  to  noise  ratio,  frequently  outperforming 
edge-based  techniques. 

For  example,  in  the  work  by  Chan  and  Vese,  a  flow  is 
proposed  which  optimally  separates  the  first  moments  of 
the  intensity  distributions  [4],  In  more  recent  work,  Rathi 
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et  al.  demonstrate  a  method  based  on  the  Bhattacharyya 
distance  for  separating  entire  distributions  [30].  In  both  of 
these  cases,  features  from  the  entire  interior  of  the  curve  are 
compared  against  features  from  the  entire  exterior. 

In  this  present  work,  we  propose  a  region-based  algo¬ 
rithm  for  the  segmentation  of  neural  fiber  bundles  from  dif¬ 
fusion  weighted  magnetic  resonance  imagery  (DW-MRI). 
Specifically,  we  describe  why  classical  approaches  (i.e. 
those  which  compare  features  across  the  full  interior  with 
features  from  the  full  exterior)  may  not  be  well-suited  for 
DW-MRI  fiber  bundle  segmentation.  Then,  we  explain  how 
one  can  leverage  the  results  of  optimal  or  geodesic  path  al¬ 
gorithms  to  locally  constrain  region-based  approaches  in 
such  a  manner  which  will  both  retain  the  beneficial  at¬ 
tributes  of  region-based  methods  while  also  handling  the 
challenges  posed  by  DW-MRI  data.  Starting  from  an  op¬ 
timal  path  (or  anchor  tract),  a  fiber  bundle  is  segmented 
using  a  Bayesian  framework.  The  priors  are  based  on 
anatomical  knowledge  of  the  bundle  being  segmented,  for 
instance,  a  simple  nonlinear  anatomically  derived  function 
of  the  distance  to  the  anchor  tract  works  well  for  the  cin¬ 
gulum  bundle.  The  likelhoods  are  based  on  local  mea¬ 
sures  of  tensor  compatibility  (local  uniformity),  adapting  a 
Chan  and  Vese  approach  to  active  contours  without  edges. 
The  Bayesian  formulation  is  cast  as  an  energy  minimization 
problem  which  is  solved  using  a  greedy  flood  fill  motivated 
algorithm. 

We  now  briefly  describe  the  remainder  of  this  paper. 
First,  in  Section  2,  we  provide  a  literature  review  and  back¬ 
ground  of  tractography  and  fiber  bundle  segmentation  algo¬ 
rithms.  Second,  in  Section  3,  we  motivate  our  interest  in 
applying  this  algorithm  to  the  segmentation  of  the  cingu¬ 
lum  bundle.  Third,  in  Section  4,  we  describe  the  algorithm 
for  locally  constraining  the  region-based  method.  Fourth, 
in  Section  5,  we  provide  initial  results  on  the  segmentation 
of  the  cingulum  bundle  using  a  simplistic  implementation. 
Finally,  in  Section  6,  we  provide  an  extensive  explanation 
of  how  these  ideas  and  results  may  be  adapted  for  use  in  a 
variety  of  implementations  and  algorithms. 
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2.  Background 

Since  the  advent  of  diffusion  weighted  magnetic  reso¬ 
nance  imaging,  a  great  deal  of  research  has  been  devoted 
to  finding  and  characterizing  neural  connections  between 
brain  structures.  Image  resolution  is  typically  high  enough 
so  that  major  white  matter  tracts,  or  bundles  of  densely 
packed  axons,  are  several  voxels  in  cross-sectional  diam¬ 
eter  [20].  The  goal  of  tractography  algorithms  is  to  segment 
these  fiber  bundles  from  the  DW-MRI  datasets. 

Early  tractography  methods  were  based  on  streamlines 
which  employed  local  decision-making  based  on  the  princi¬ 
pal  eigenvector  of  diffusion  tensors  [19,  33,  2,  5]1. 

In  these  techniques,  tracts  are  propagated  from  a  start¬ 
ing  point  until  the  tracts  reach  some  termination  criterion. 
Due  to  the  local  decision-making  process,  these  methods 
have  been  shown  to  perform  poorly  in  noise  and  often  stop 
prematurely.  These  techniques  do  not  provide  a  measure  of 
connectivity  for  the  resulting  tracts.  Furthermore,  several  of 
these  methods  do  not  use  the  full  tensor,  reducing  the  data 
to  the  principal  eigenvectors,  and  subsequently  are  unable 
to  handle  fiber  crossings,  branchings,“kissings,”  etc. 

Despite  the  shortcomings  of  this  approach,  due  to  its 
ease-of-use,  streamlining  has  quickly  become  the  most  pop¬ 
ular  method  for  fiber  segmentation.  To  infer  fiber  bundles 
from  streamline  tractography  results,  several  groups  have 
successfully  worked  on  methods  for  fiber  clustering.  The 
goal  of  clustering  is  to  capture  group  behavior  of  a  popula¬ 
tion  of  streamlines  and  to  use  this  group  behavior  to  drive 
fiber  bundle  segmentation.  The  end  result  of  clustering  al¬ 
gorithms  has  been  shown  to  accurately  capture  many  neural 
fiber  bundles,  see  for  example  [22,  18], 

Recently,  another  line  of  work  has  emerged  which  seeks 
to  avoid  the  use  of  the  problematic  streamlines.  Tractogra¬ 
phy  advances  have  been  made  which  provide  full  brain  opti¬ 
mal  connectivity  maps  from  predefined  seed  regions.  These 
methods  are  more  robust  to  noise  and  depending  upon  the 
underlying  metric,  may  be  able  to  more  fully  use  the  com¬ 
plete  DW-MRI  data.  These  approaches  can  be  subdivided 
into  stochastic  and  energy-minimization  approaches. 

Stochastic  approaches  produce  probability  maps  of  con¬ 
nectivity  between  a  seed  region  and  the  rest  of  the  brain. 
Parker  et  al.  developed  PICo,  a  probabilistic  index  for  stan¬ 
dard  streamline  techniques  [24],  Perrin  et  al.  presented 
probabilistic  techniques  for  untangling  fiber  crossings  using 
q-ball  fields  [26].  In  other  work,  Friman  et  al.  proposed  a 
method  for  probabilistically  growing  fibers  in  a  large  num¬ 
ber  of  random  directions  and  inferring  connectivity  from  the 
resulting  percentages  of  connections  between  seed  and  tar- 

1  The  diffusion  tensor  is  one  of  the  simplest  diffusion  models.  It  is 
estimated  from  a  set  of  diffusion  weighted  images,  each  probing  the  water 
diffusion  in  a  different  spatial  direction.  In  the  three-dimensional  case  the 
diffusion  tensor  is  a  3  X  3  symmetric,  positive  definite  tensor.  For  details 
see  [3]. 


get  regions  [7].  While  providing  a  measure  of  connectivity 
between  brain  regions,  these  stochastic  approaches  do  not 
provide  an  explicit  segmentation  of  the  fiber  bundle  itself 
and  often  do  not  explicitly  provide  the  optimal  connection 
between  regions  of  the  brain. 

Energy-minimization  techniques  have  also  been  devel¬ 
oped.  Parker  et  al.  proposed  fast  marching  tractogra¬ 
phy  which  minimizes  an  energy  based  on  both  the  posi¬ 
tion  and  direction  of  the  normal  to  a  propagating  front  [25]. 
O’Donnell  et  al.  cast  the  tractography  problem  in  a  geo¬ 
metric  framework  finding  geodesics  on  a  Riemannian  man¬ 
ifold  based  on  diffusion  tensors  [23] .  Similarly,  Prados  et  al. 
and  Fenglet  et  al.  demonstrated  a  Riemannian  based  tech¬ 
nique,  GCM  (Geodesic  Connectivity  Mapping),  for  com¬ 
puting  geodesics  using  a  variant  of  fast  marching  methods 
adapted  for  directional  flows  [29,  15].  Jackowski  et  al. 
find  Riemannian  geodesics  using  Fast  Sweeping  methods 
as  given  by  Kao  et  al.  [8,  11,  12].  Pichon  et  al.  and  Mel- 
onakos  et  al.  use  the  more  general  Finsler  metric  to  find 
optimal  connections  [27,  28,  17,  16].  Finally,  Fletcher  et  al. 
propose  a  new  Hamilton-Jacobi-Bellman  numeric  solver  on 
the  graphics  processing  unit  to  find  Riemannian  geodesics 
in  near  real-time  speeds  [6],  In  each  of  these  cases,  an  opti¬ 
mal  path  is  found  which  represents  the  best  connection  be¬ 
tween  the  two  regions  under  the  given  metric. 

3.  The  Cingulum  Bundle 

In  this  section,  we  motivate  the  problem  of  segment¬ 
ing  the  cingulum  bundle.  The  cingulum  bundle  is  a  5-7 
mm  in  diameter  fiber  bundle  that  interconnects  all  parts  of 
the  limbic  system.  It  originates  within  the  white  matter  of 
the  temporal  pole,  and  runs  posterior  and  superior  into  the 
parietal  lobe,  then  turns,  forming  a  ”ring-like  belt”  around 
the  corpus  callosum,  into  the  frontal  lobe,  terminating  ante¬ 
rior  and  inferior  to  the  genu  of  the  corpus  callosum  in  the 
orbital-frontal  cortex  [32].  Moreover,  the  cingulum  bun¬ 
dle  consists  of  long,  association  fibers  that  directly  connect 
temporal  and  frontal  lobes,  as  well  as  shorter  fibers  radiat¬ 
ing  into  their  own  gyri.  The  cingulum  bundle  also  includes 
most  afferent  and  efferent  cortical  connections  of  cingulate 
cortex,  including  those  of  prefrontal,  parietal  and  temporal 
areas,  and  the  thalamostriatae  bundle.  In  addition,  lesion 
studies  document  a  variety  of  neurobehavioral  deficits  re¬ 
sulting  from  a  lesion  located  in  this  area,  including  akinetic 
mutism,  apathy,  transient  motor  aphasia,  emotional  distur¬ 
bances,  attentional  deficits,  motor  activation,  and  memory 
deficits.  Because  of  its  involvement  in  executive  control  and 
emotional  processing,  the  cingulum  bundle  has  been  inves¬ 
tigated  in  several  clinical  populations,  including  depression 
and  schizophrenia.  Previous  studies,  using  diffusion  tensor 
imagery,  in  schizophrenia,  demonstrated  decrease  of  frac¬ 
tional  anisotropy  in  the  anterior  part  of  the  cingulum  bun¬ 
dle  [13,  34],  at  the  same  time  pointing  to  the  technical  lim- 


itations  restricting  these  investigations  from  following  the 
entire  fiber  tract. 

4.  The  Algorithm 

In  this  section,  we  present  our  method  for  applying  local 
constraints  to  region-based  flows. 

4.1.  Motivation  for  Local  Constraints 

An  implicit  assumption  of  classical  (i.e.,  those  which 
compare  features  across  the  full  interior  with  features  from 
the  full  exterior)  region-based  approaches  is  that  the  en¬ 
tire  interior  of  the  contour  contains  fairly  homogeneous  fea¬ 
tures,  such  as  mean  intensity.  Under  this  assumption,  these 
algorithms  proceed  by  evolving  the  closed  curve  or  surface 
to  minimize  an  energy  defined  over  these  features. 

However,  if  there  are  no  homogeneous  features  across 
the  entire  interior  or  exterior  of  the  object  of  interest,  it 
becomes  difficult  to  define  a  region-based  approach  which 
will  accurately  segment  the  image.  For  instance,  in  the  case 
of  the  cingulum  bundle  which  curves  around  the  ventricles, 
the  tensors  across  the  fiber  bundle  vary  in  both  anisotropy 
and  orientation  across  the  length  of  the  bundle,  as  shown 
in  Figure  1 .  In  this  sagittal  view,  we  see  that  it  is  difficult 
to  define  a  feature  on  the  space  of  tensors  which  uniquely 
separates  the  entire  interior  of  the  cingulum  bundle  from  the 
exterior.  However,  we  also  notice  that  the  tensor  shape  and 
anisotropy  vary  smoothly  across  the  bundle.  Hence,  locally 
across  the  fiber  one  can  define  tensor  features  which  are  dis¬ 
tinguishable  from  the  exterior. 

4.2.  Prior  Work 

Surface  evolution  approaches  have  been  described  for 
fiber  bundle  segmentation.  Rousson  et  al.  [31]  use  a  multi¬ 
variate  Gaussian  distribution  of  the  tensor  components  in  a 
geodesic  active  region  model  to  drive  a  surface  evolution 
towards  the  segmentation  of  fiber  bundles.  The  method  is 
applied  to  the  segmentation  of  the  corpus  callosum,  but  is 
unable  to  fully  capture  its  curved  character  as  discussed  by 
the  authors.  In  a  follow-up  paper  [14]  a  similar  segmen¬ 
tation  framework  in  combination  with  a  geodesic  distance 
between  tensors  is  shown  to  yield  superior  segmentation  re¬ 
sults,  in  particular,  when  segmenting  curved  fiber  bundles. 
Jonasson  et  al.  propose  two  different  ways  to  address  the 
segmentation  of  curved  fiber  bundles  in  a  surface  evolution 
setting:  (i)  a  local  approach  [9],  where  the  surface  evolution 
speed  is  influenced  by  the  similarity  of  a  tensor  in  com¬ 
parison  to  its  interior  neighbors,  and  (ii)  a  region-based  ap¬ 
proach,  where  the  similarity  measure  is  based  on  the  no¬ 
tion  of  a  most  representative  tensor  within  the  segmented 
region  [10].  In  the  latter  case,  capturing  highly  curved  fiber 
bundles  will  be  problematic.  In  both  cases  the  segmenta¬ 
tion  algorithm  is  combined  with  a  surface  regularization  to 


prevent  leaking.  The  approach  proposed  in  this  paper  is 
related  to  Jonasson’s  work  [9]  in  as  much  as  it  uses  local 
tensor  similarities  to  drive  the  segmentation,  however,  no 
surface  evolution  is  used  and  a  tensor  similarity  measure 
is  combined  with  prior  information  as  given  by  an  initially 
computed  anchor  tract  (also  preventing  large-scale  leaking). 
The  extension  of  the  approach  proposed  in  this  paper  (see 
Section  6)  can  be  seen  as  complementary  to  the  method  by 
Lenglet  et  al.  [14].  Instead  of  disentangling  tensor  shape 
and  orientation  through  an  appropriate  tensor  distance  (and 
statistic)  the  anchor  tract  may  be  used  to  warp  the  space 
initially,  thus  effectively  removing  large  orientation  differ¬ 
ences2.  Further,  due  to  the  absence  of  a  surface  evolution, 
our  approach  is  computationally  very  efficient. 

4.3.  Bayesian  Framework 

In  this  section,  we  describe  how  the  algorithm  can  be  for¬ 
mulated  in  a  Bayesian  framework.  We  follow  the  approach 
by  Mumford  [21]  and  cast  the  Bayesian  estimation  problem 
into  an  energy  minimization.  The  probability  of  observing 
the  classification  C ,  consisting  of  points  belonging  to  the 
fiber  and  points  belonging  to  the  background  given  the  ten¬ 
sor  information  T  is  (using  Bayes’  formula): 

P(C\T)  =  P{T^C)  ~  P(T\C)P(C),  (1) 

where  C  is  an  element  of  the  set  of  all  possible  assign¬ 
ments  of  voxels  to  the  fiber  and  the  background  respectively, 
p(T\C)  is  the  likelihood  of  observing  T  given  the  classifi¬ 
cation  C  and  p(C)  is  the  prior.  By  taking  the  logarithm 
on  both  sides  and  noting  that  p(T)  is  independent  of  the 
classification  C ,  Equation  (1)  can  be  written  as  an  energy 
minimization  problem  [21] 

E(C)  =  -log  (p(C,T)) 

-  -log(p(T|C))-log(P(C)) 

=  Ed(T,C)  +  Ep(C),  (2) 

where  Ed{T,  C)  denotes  the  data  energy  and  Ep(C )  the 
prior  (or  regularization)  energy.  Instead  of  solving  the 
Bayesian  estimation  problem  (1)  directly  we  may  thus  in¬ 
stead  minimize  the  energy  (2).  Which  leaves  us  with  defin¬ 
ing  these  energies.  We  use  a  flood-fill  algorithm  approach 
that  solves  the  energy  minimization  problem  (2)  for  an  in¬ 
dividual  point  only  considering  its  local  neighborhood  N. 
In  what  follows  we  first  describe  the  continuous  setting,  to 
make  connections  with  existing  approaches,  and  then  de¬ 
scribe  the  discrete  implementation  in  the  context  of  the  pro¬ 
posed  Bayesian  flood-fill  algorithm.  Given  the  local  neigh¬ 
borhood  N  of  a  point  x  we  want  to  decompose  it  into  a  sub- 
region  belonging  to  the  fiber  and  a  subregion  belonging  to 

2  Our  approach  may  also  be  combined  with  the  method  proposed 
in  [14]. 


Figure  1.  Example  of  the  need  for  local  constraints  on  region-based  segmentation  algorithms  which  attempt  to  segment  the  cingulum 
bundle.  Notice  that  tensor  anisotropy  and  orientation  vary  across  the  length  of  the  cingulum  bundle. 


the  background.  The  goal  of  our  algorithm  is  to  make  each 
of  these  two  subregions  individually  as  uniform  as  possi¬ 
ble,  while  at  the  same  time  using  anatomically  meaningful 
prior  information.  The  prior  information  is  encoded  based 
on  the  distance  of  the  pre-computed  anchor  tract,  which  is 
the  lowest  cost  path  connecting  two  maximally  spaced-out, 
pre-defined  regions  of  interest  of  the  fiber  bundle  of  interest 
(in  our  case  the  cingulum  bundle).  Specifically,  we  choose 
p(C)  as 

p(x)  =  sG(d(x)),  (3) 

where  d(x)  is  the  distance  of  point  x  from  the  anchor  tract 
and 

!1  for  \r\  >  Umin 

5  for  Umin  <  \r  I  <  Umax 

0  otherwise, 

where  Ga  is  a  Gaussian  with  standard  deviation  a  and  * 
is  the  convolution  operator;  pmin  and  p,max  are  set  to  the 
range  of  expected  radius  values.  Note,  that  the  prior  could 
also  be  replaced  by  a  probabilistic  atlas.  Equation  (3)  de¬ 
scribes  an  initial  zone  of  high  fiber  confidence  close  to  the 
anchor  tract,  a  transitioning  region  (where  p(C)  =  1/2) 
where  the  prior  information  will  not  be  used3,  and  an 
anatomically  implausible  region,  where  the  prior  probabil¬ 
ity  decreases  to  zero.  The  prior  energy  is  then  defined  as 

EP(C)  =  1  -  p(x)  dfl  +  J  p(x)  <m^J  ,  (4) 

where  Nf  is  the  region  belonging  to  the  fiber  Nb  is  the  re¬ 
gion  belonging  to  the  background  and  |  •  |  denotes  cardinal¬ 
ity,  i.e.,  \N\  is  the  volume  of  the  neighborhood.  Given  a 

’If  p(C)  =  1/2,  the  prior  energy  (4)  is  independent  of  assigning  the 
candidate  flood-fill  point  to  the  fiber  or  the  background. 


measure  of  uniformity  D  :  T  x  Sn  i— >  Rj  mapping  from 
the  space  of  tensors  T  and  the  space  of  neighborhood  sets 
of  tensors  Sn  3  T(N)  :=  {T(x)  £  T|x  £  N}  to  a  non¬ 
negative  real  value,  we  write  the  data  energy  as 

Ed{T,c)=w\{L  DvrW'T(Nf»dn 

+  [  D(T(x),T(Nb))  dfi ),  (5) 

JNb  J 

where  T(x)  denotes  a  tensor  at  position  x  and  T(N)  de¬ 
notes  the  set  of  tensors  in  the  region  N.  This  is  an  energy 
similar  to  the  one  proposed  by  Chan  and  Vese  [4]  for  the 
segmentation  of  intensity  images4.  Note,  however,  that  in¬ 
stead  of  using  this  energy  globally  to  perform  tensor  seg¬ 
mentation,  we  are  proposing  to  use  this  energy  in  a  local 
neighborhood  to  make  a  local  decision  for  a  flood-fill  al¬ 
gorithm,  thus  avoiding  global  tensor  orientation  issues  for 
strongly  curving  fiber  bundles.  To  minimize  this  energy  in 
the  discrete  flood-fill  setting,  we  simply  compute  the  dif¬ 
ference  of  the  energies  when  adding  the  voxel  in  question 
to  either  the  fiber  (resulting  in  energy  E f)  or  to  the  back¬ 
ground  (resulting  in  energy  Eb).  The  difference  of  these 
energies  A E  =  Eb  —  E*  corresponds  to  a  discretized  gra¬ 
dient.  Since  our  goal  is  to  minimize  the  overall  energy,  a 
voxel  x  will  be  added  to  the  set  of  fiber  voxels  if  A E  >  0. 
All  integrals  in  Equations  (4)  and  (5)  become  sums  in  the 
discretization.  Many  uniformity  measures  are  possible  (see 
for  example  [10,  9,  1]  for  some  ideas  on  how-to  compare 
tensors),  we  constructed  a  simple  one  based  on  fractional 


4To  favor  “smooth”  discrete  boundaries,  a  local  boundary  length  term 
can  be  added. 


anisotropy  and  the  major  diffusion  direction: 

£(T(x),T(iV))  =  j(%(T(x),rffl) 

+  Dei(T(x),T(iV))), 


(iv)  Repeat  from  step  (ii)  until  no  more  new  fiber  voxels 
are  found. 

Using  the  Bayesian  framework,  the  outward  propagating 
front  stops  once  the  Bayesian  detection  threshold  is  reached, 
i.e.,  once  all  boundary  voxels  are  in  locally  minimal  energy 
configurations. 


where 

Dfa(T(x),T(N))  =  | FA(T(x))  -  FA(T(N))\ 
measures  the  uniformity  in  fractional  anisotropy  and 


£>ei(T(x),T(JV))  =  1  -  y/FA(T(x))FA(T(N)) 


x  ei (T(x)) 


t  A1(T(iV))e1(r(7V)) 
||A1(T(iV))e1(T(7V))|| 


measures  the  uniformity  in  direction.  Fractional  anisotropy 
(FA)  is  defined  as  [3] 


pA=  /3  -  A2)2  +  (A!  -  A3)2  +  (AT^AsF 

V  2  y/A?  +  A|  +  A| 

where  ei(T)  denotes  the  major  unit  eigenvector  of  the  ten¬ 
sor  T,  A i(T)  its  eigenvalues  (with  Ai  >  A2  >  A3  >  0), 
and  the  overhead  bar  signifies  the  mean5.  Dei  is  scaled  by 
fractional  anisotropy  to  discard  tensors  that  are  close  to  be¬ 
ing  isotropic,  since  in  these  cases  eigenvector  computations 
become  numerically  problematic.  The  continuous  approach 
could  alternatively  be  implemented  using  fast  marching  or 
level  sets.  In  this  work,  we  use  a  very  simple  flood-fill  ap¬ 
proach  which  propagates  away  from  the  anchor  tract.  Cer¬ 
tainly  other  methods  would  offer  a  more  continuous  and  nu¬ 
merically  accurate  approach.  However,  our  simple  flood-fill 
implementation  is  sufficient  as  a  proof-of-concept. 

The  algorithm  proceeds  in  the  following  steps: 

(i)  Declare  all  voxels  on  the  anchor  tract  as  fiber  voxels. 

(ii)  Consider  all  6-connected  neighbors  to  the  fiber  voxels 
that  are  not  fiber  voxels  themselves  as  candidate  vox¬ 
els. 

(iii)  Decide  whether  a  candidate  voxel  should  belong  to  the 
fiber  based  on  the  simple  local  energy  minimization 
described  above  (where  the  neighborhoods  Nf  and  Nt, 
are  given  by  the  voxels  in  the  current  neighborhood 
N  that  already  belong  to  the  fiber  or  are  so  far  classi¬ 
fied  as  background  respectively).  If  a  candidate  voxel 
should  be  part  of  the  fiber  according  to  the  local  energy 
minimization,  add  it  as  a  fiber  voxel. 

5FA  may  be  computed  directly  from  the  tensor  components  without 
computing  the  tensor  eigenvalues  first  [3]. 


5.  Experiments 

In  this  section,  we  show  results  of  the  algorithm  applied 
to  DW-MRI  datasets  of  51  sampling  directions.  We  used  the 
Finsler  tractography  method  proposed  by  Melonakos  et  al. 
to  compute  the  anchor  tracts  or  optimal  paths  between  two 
input  seed  regions  [16].  The  seed  regions  were  manually 
segmented,  one  under  the  anterior  tip  of  the  ventricles  and 
the  other  under  the  posterior  tip  of  the  ventricles. 

Using  precomputed  anchor  tracts,  we  were  able  to  con¬ 
struct  our  priors  using  the  function  shown  in  Figure  2,  as 
previously  described  (mean  radius  r  =  3  mm,  /jm,n  = 
\r  =  1.5  mm,  nmax  =  § r  =  4.5  mm,  a  =  |r  =  |).  Ap¬ 
plying  this  function  to  a  distance  map  from  the  anchor  tract, 
the  prior  image  is  as  shown  on  the  left  side  of  Figure  3.  The 
white  colored  area  is  where  the  uniform  priors  are  centered 
on  the  mean  value  of  the  cingulum  bundle  radius,  which  we 
take  to  be  3  mm  as  described  in  Section  3.  In  the  middle  of 
Figure  3,  we  show  the  likelihood  energy  gradient  computed 
from  the  evolution  (where  positive  values  are  likely  to  be¬ 
long  to  the  bundle  and  negative  values  are  not  likely  to  be¬ 
long  to  the  bundle).  Notice  how  the  likelihood  energy  func¬ 
tion  captures  an  appropriate  boundary  across  a  majority  of 
the  cingulum  fiber  bundle.  The  orientation  dependent  terms 
had  the  strongest  influence  on  the  inferior  edge  against  the 
corpus  callosum.  The  anisotropy  dependent  terms  had  the 
strongest  influence  on  the  superior  edge.  On  the  right  side 
of  Figure  3,  we  show  the  posterior  energy  gradient,  which 
results  from  the  combination  of  likelihood  energy  and  prior 
energy  terms. 

In  Figure  4,  we  show  a  3D  model  view  of  the  result¬ 
ing  segmentation.  Then,  in  Figure  5,  we  show  three  sepa¬ 
rate  time  steps  in  the  flood-fill  evolution.  The  first  column 
is  at  1  iteration,  the  second  column  is  at  3  iterations,  and 
the  final  column  is  at  18  iterations-where  all  three  methods 
had  converged.  The  top  row  shows  the  evolution  using  only 
the  priors.  Notice  how  the  result  is  a  smooth  tube  exactly 
matching  the  prior  that  is  too  wide  for  this  individual  and 
ends  up  overlapping  proximal  anatomy,  such  as  the  corpus 
callosum.  The  middle  row  shows  the  evolution  using  only 
the  likelihoods.  While  this  result  appropriately  captures  the 
majority  of  the  bundle,  it  is  subject  to  a  few  leaks  as  shown. 
The  bottom  shows  the  evolution  using  the  Bayesian  combi¬ 
nation  of  the  likelihoods  and  priors.  This  result  shows  an 
appropriate  combination  of  the  likelihood  boundary  stop¬ 
ping  and  the  prior  leakage  constraints. 


Figure  3.  The  prior  energy  (left),  likelihood  energy  (middle),  and  posterior  energy  (right). 


Figure  4.  A  3D  view  of  the  result. 


Figure  5.  Front  evolution  time  steps:  Top  Row  is  the  evolution  with  only  the  priors.  Middle  Row  is  the  evolution  with  only  the  likelihoods. 
Bottom  Row  is  the  evolution  from  the  Bayesian  inclusion  of  both  the  likelihoods  and  priors. 


Figure  2.  The  prior  profile:  Blue  is  the  initial  step  function.  Red  is 
the  actual  profile  after  smoothing.  Note  the  region  of  uniform  pri¬ 
ors  (0.5),  centered  around  the  clinically  defined  mean  fiber  radius. 


We  also  note  that  the  few  parameters  used  in  this  method 
can  be  chosen  given  anatomical  information  about  the  mean 
radius  of  the  fiber  bundle.  The  prior  energy  function  is  only 
dependent  upon  this  parameter,  as  mentioned  previously. 
Also,  the  neighborhood  size  is  chosen  to  be  large  enough 
so  that  at  least  20%  of  the  neighborhoods  on  the  first  it¬ 
eration  include  voxels  exterior  to  the  fiber  bundle.  In  this 
case,  we  chose  a  neighborhood  radius  of  7  mm.  No  other 
parameters  were  needed  in  this  computation. 

6.  Future  Work  and  Conclusions 

This  paper  proposed  a  novel  segmentation  method  for 
diffusion  tensor  images.  The  approach  is  based  on  a 
Bayesian  region  growing,  where  the  prior  depends  on  the 
distance  to  a  pre-computed  anchor  tract.  The  anchor  tract 
is  given  by  the  optimal  path  in  a  Finsler  metric  (though  any 
other  robust  method  giving  a  representative  fiber  path  could 
be  used),  utilizing  the  full  diffusion  profile.  The  likelihood 
is  determined  based  on  the  consistency  of  a  candidate  voxel 
with  its  neighbors  that  are  already  part  of  the  segmentation, 
(i.e.  the  likelihood  is  dynamically  updated  as  the  region 
is  growing.)  The  Bayesian  combination  of  likelihood  and 
prior  allows  for  a  balanced  combination  of  local  consistency 
and  distance  from  the  optimal  path,  which  also  inhibits  seg¬ 
mentation  leakage.  The  approach  is  computationally  effi¬ 
cient. 

Region-based  segmentation  algorithms  have  been  highly 
successful  in  segmenting  uniform  (in  a  given  measure)  re¬ 
gions.  Translating  this  global  region-based  approach  to  dif¬ 
fusion  weighted  imaging  for  the  segmentation  of  fiber  bun¬ 
dles  is  challenging,  since  it  is  not  obvious  how  to  define  a 
criterion  to  incorporate  shape  and  directional  information 
over  an  entire  region.  In  particular,  many  fiber  bundles  in 
the  brain  curve  strongly  (e.g.  the  cingulum  bundle,  the  arcu¬ 
ate  fasciculus,  the  corpus  callosum).  Figure  6  shows  some 


Figure  6.  Tensors  aligned  along  the  anchor  tract  (top)  and  tensors 
aligned  along  the  anchor  tract  warped  to  a  straight  line  (bottom). 
Warping  the  tensors  based  on  the  geometry  of  the  anchor  tract 
greatly  simplifies  the  tensor  segmentation  problem. 


exemplary  tensors  aligned  along  an  anchor  tract.  Warping 
the  anchor  tract  to  a  straight  line  may  greatly  simplify  the 
design  of  region-based  statistics,  by  “flattening”  the  geom¬ 
etry  to  remove  large-scale  deformation.  Further,  by  estab¬ 
lishing  correspondences  between  the  anchor  tract  and  the 
rest  of  the  domain,  other  interesting  neighborhoods  may  be 
defined.  This  will  be  the  topic  of  future  research. 

Many  more  extensions  are  conceivable,  such  as  the  use 
of  more  sophisticated  distance  and  similarity  measures. 
The  locally-constrained  method  proposed  as  well  as  global 
region-based  segmentation  methods  will  benefit  from  sim¬ 
ilarity  metrics  using  the  complete  tensor  information.  In 
particular,  more  suitable  tensor-based  statistics  may  be  ex¬ 
plored  in  this  framework,  such  as  those  shown  by  Lenglet 
et  al.  [14].  Further,  a  continuous  formulation,  based  on  a 
variant  of  the  Eikonal  equation  or  as  a  complete  surface 
evolution  (with  the  easy  possibility  of  directly  integrating 
smoothing  terms)  will  be  desirable,  and  validations  with  re¬ 
spect  to  manual  segmentations  should  be  performed. 
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